Radiolysis‐Driven Evolution of Gold Nanostructures – Model Verification by Scale Bridging In Situ Liquid‐Phase Transmission Electron Microscopy and X‐Ray Diffraction

Abstract Utilizing ionizing radiation for in situ studies in liquid media enables unique insights into nanostructure formation dynamics. As radiolysis interferes with observations, kinetic simulations are employed to understand and exploit beam‐liquid interactions. By introducing an intuitive tool to simulate arbitrary kinetic models for radiation chemistry, it is demonstrated that these models provide a holistic understanding of reaction mechanisms. This is shown for irradiated HAuCl4 solutions allowing for quantitative prediction and tailoring of redox processes in liquid‐phase transmission electron microscopy (LP‐TEM). Moreover, it is demonstrated that kinetic modeling of radiation chemistry is applicable to investigations utilizing X‐rays such as X‐ray diffraction (XRD). This emphasizes that beam‐sample interactions must be considered during XRD in liquid media and shows that reaction kinetics do not provide a threshold dose rate for gold nucleation relevant to LP‐TEM and XRD. Furthermore, it is unveiled that oxidative etching of gold nanoparticles depends on both, precursor concentration, and dose rate. This dependency is exploited to probe the electron beam‐induced shift in Gibbs free energy landscape by analyzing critical radii of gold nanoparticles.

electron beam containing an isotropic distribution of all contributing species is simulated.
Under homogeneous irradiation, a quasi-closed system can form within the irradiated volume with diffusion of chemical species towards non-irradiated volumes being negligible. [32,47] In this voxel, the chemical interplay is described by a set of coupled ordinary differential equations (ODEs). In the kinetic model, each ODE depicts the changes in concentration c over time t of a reactant i [21] : In equation (S1), the first term relates to the physical interaction between the ionizing radiation and the liquid phase, with liquid density ρ, Faraday constant F, dose rate ψ, and generation value Gi, which can be approximated to zero for solvated species at small concentrations. [37] The generation values used in this work are listed in Table S1. The second and third terms account for creation/consumption of species i during a chemical reaction with the rate constant k.  [21] X-Ray irradiation / (Molecules/100 eV) [74] -5.68 -4.64 In the past, such simulations have been achieved by manually creating the set of ordinary differential equations to solve equation (S1). This is a cumbersome and highly complicated process, prone to hard-to-detect errors. Consequently, extensions of the reaction set comprising pure water are sparse and (mostly) limited to a few additional reactions.

H2O
A routine solely relying on open-source software is developed to overcome these limitations.
The script based on Python 3.7, NumPy [48] , Matplotlib [49] , pandas [50] , and SciPy [51] takes a plaintext file of the chemical reaction set as input and automatically creates the corresponding matrix of coupled differential equations. Its workflow is sketched in Figure 1(a) of the main manuscript.
As illustrated in Figure 1(b), this is achieved by using natural language processing via regular expressions (regex). Furthermore, the software is equipped with a tool to check the input reactions for consistent charge and atom balance. Next, the ODE matrix is translated into an array to provide suitable input for the solve_ivp function of SciPy, where each row follows the syntax of equation (1). This is illustrated in Figure 1(c).
By using the Fortran-based "LSODA" solver of SciPy, the solution of the coupled ODE-system is determined numerically. An exemplary output is shown Figure 1(d) for the model described below, consisting of 42 different chemical species distributed over 184 reactions.
The electron flux density ϕ is used to account for experimentally accessible input variables. Via ϕ = n e / ( ), with the elementary charge e, the number of electrons n per time interval t (here: 1 s), and area (here: 1 nm²), it is converted into the dose rate experienced during the experiment: [18] = (1 + ) with d as liquid thickness (d ≤ 100 nm in this study), λ as inelastic mean free path (380 nm at 300 kV acceleration voltage [52] ), and S as stopping power (2.36 MeV(cm)²g -1 at 300 keV electrons in water [53] ).

S1.2 Validation
The tool has been validated against the MATLAB implementation of the reaction set of Schneider et al., [21] and shows excellent agreement in both temporal ( Figure S1(a)) and steadystate evolution (Figure S1(b)).
The code is available at GitHub at https://github.com/BirkFritsch. It is held in a modular fashion to simplify possible extensions and incorporations into third-party software. The code can be easily executed using Anaconda, or any other Python environment with the required third-party modules installed. A list of the dependencies is given on GitHub.
Furthermore, it is extended by a graphical user interface (GUI) allowing for simulation of reaction kinetics without any programming knowledge ( Figure S2). The GUI is also available on GitHub. The tool incorporates an option to check for errors in the reaction text file by verification of mass-and charge balance for each individual reaction and a scan for doublings.
Doublings are reactions where an identical set of educts yields an identical set of products.
Furthermore, the tool offers the option to verify that mass was conserved during the kinetic simulation. This differs significantly from the mass balance check of the reaction text file, because here, also the generation value balance is crucial. Figure S6 provides an exemplary output for the reaction set described below.
Figure S1: Comparison of water radiolysis simulation using the tool presented here and the MATLAB-script published by Schneider et al. [21] . (a) Temporal simulation at 20 MGy s -1 showing only the radiolysis products. (b) Steady-state concentrations including the solvent concentration. Figure S2: Screenshot of the input text files, the GUI window, and its main functionalities.
The reaction set consists of three columns separated by a semicolon ";" (see Figure S2 for an example of first eleven reactions used in this work): • The left-hand side is the most important part, as it describes the chemical reaction: o Educts and products are separated by " -->". o Arbitrary names for agents can be used, if they are used consistently throughout the whole set. This in turn also means that every name must be unique. Note that H* and H would be two different species! Extra attention has to be paid when copying content from another source, e.g. a pdf file, etc.. Only because symbols look similar, they are not necessarily the same ASCII symbol. A major issue is "-" and "-". ▪ WARNING: While the kinetic simulation itself does not require meaningful names; the error handling (partially) requires a chemistry-like notation. Hence, we suggest using the following convention that leans on the usual nomenclature for chemical formulae: • Every chemical element starts with a capital letter, e.g. H, Au.
• Indices denoting the amount of an element in a chemical formula are written as Arabic numbers in-line. (e.g. H2O, Cl2, Au), whereas an index of "1" can be omitted. • Charges are written at the end of a sum formula. In case of a single charge, a sole "+" or "-" suffices. If more charges are present, this should be indicated by a "^" followed by an Arabic number (e.g. Cl-, H+, Au2Cl6^2-for Cl -, H + , or Au2Cl6 2-). • Parentheses indicating oxidation states are fine (e.g. Au(III)Cl4-). o Stoichiometric factors are written in front of the reactant separated by a space " ". E.g. 2 H2O. o Agents and their stoichiometric factors are separated by " + ". Take care to put the space in front and after the "+" so that the code recognizes it as an operator. • The right-hand side is providing the reaction rate constants in mol -n+1 dm 3(n-1) s -1 with n as reaction order. You can abbreviate large numbers by using scientific notation e.g. 5e9 (= 5 ⋅ 10 9 ) instead of 5000000000, or 1.3e-1 instead of 0.13. o It is important, however, to use "." as decimal spacer and do not type any spacer symbols such as " " or ",". o Do not type a unit here, just provide a number. • These would be valid lines in the code: o 2 eh-+ 2 H2O --> H2 + 2 OH-; 5.5e9 ; a comment. o H+ + eh---> H ; 0 ; just another comment. • Do not leave blank lines in between the reaction set.
• Do not write headings in the reaction set.
The second file is related to the initial concentrations and G-values of the species if known (see again Figure S2 for an example of first four radiolytic species used in this work): • Here, you can enter a starting concentration in M and a constant G-value if required.
• It is not necessary to list all species here. If it is not listed, the program displays a warning and automatically sets both, G value and the initial concentration of this species to 0. Again, pay attention to use the identical syntax for your agents in the file as in the reactions. • The file consists of three columns separated by a tab stop.
o The first one lists the species name identical to the used one in the reaction set file, the second one the initial concentration in M, and the third one the G-value in molecules/100eV. • Do not comment in this file.
• Do not leave blank lines in between. Store the reaction set as ".txt" file, e.g. "concentrations.txt".
When creating a reaction set it is important to ensure that every species included exhibits at least one formation and one decay pathway. If the former is missing, it will not be formed. In the latter case, its concentration will inevitably rise until one precursor species is exhausted.
However, as redox reactions are in principle reversible such a simulation is usually meaningless if the counter reaction is not regarded.
In our experience, it is possible to set up a simulation for a novel reaction set within short notice, if all required reactions and parameters are known. As elucidated by Schneider et al., [21] the NIST solution kinetics database (https://kinetics.nist.gov/solution/) provides a good starting point to survey possible radiation-chemistry reactions.
A more comprehensive manual of the GUI and working examples including the HAuCl4 set used in this work (see Table S2 below) can be found on GitHub.

S1.3 Kinetic modeling of HAuCl4 solutions
During this work, a comprehensive reaction set for HAuCl4 irradiation consisting of 184 reactions and 42 reactive species has been created, based on the set proposed by Schneider et al. [21] . On top, aqueous species are appended by the O-radical and O4 to account for their importance in reactions with chlorine species and the spontaneous decomposition of H2O2 and O3. The interaction of gold with radiolytic species is based on Dey et al. [54] and Ambrožič et al. [38] . The chlorine reactions mainly stem from combining the works of Kläning, Sehested, and Holcman [56] , Kelm and Bohnert [57] , and Levanov and Iskaikina [58] , while the amount of chloride ligands in the gold species proposed Dey et al. [54] , and Buxton and Sellers [59] is following Gachard et al. [60] . In accordance with Park et al. [37] , HAuCl4 was assumed to be fully deprotonated prior to irradiation. As shown in Figure S3, steady states are formed within a fraction of a second at dose rates relevant to TEM.
A graph representation of the reaction set is shown in Figure S4. In graph theory, every reactant can be represented by a node. The links (arrows) connecting the nodes represent every chemical reaction where the reactant is involved. An arrow spanning from reactant A to reactant B is to be read as 'A is an educt forming B'. Thus, for a schematic reaction A + B → C + D, four arrows would be drawn, connecting A with C, A with D, B with C, and B with D. This concept is illustrated in the legend of Figure S4. Additionally, the weight of every link represents the decadic logarithm of the respective kinetic constant times the stoichiometric factor of the educt node, contributing to the width of the link, as suggested elsewhere. [68] Figure S3: Temporal evolution of the concentration of (a) water-based and (b) gold/chlorinecontaining species of irradiation of 20 mм HAuCl4 at 1000 e -(nm²s) -1 . Figure S4: Graph representation of the utilized reaction set.
To account for the relevance of the respective reactant within the reaction set, the size of a node is representing its betweenness centrality, a measure corresponding to the direct connections of different species traversing this node. The color of a link matches the color of the node where it originates from. Graphs are created via NetworkX. [67] To account for the generation of primary species, they are accounted for by pseudo-reactions with a pseudo-reaction constant of Gi•10 12 s -1 .
To analyze the relevance of a species within the reaction set, the harmonic centrality H is used.
It can be understood as the amount of information traversing a node within a network with weighted links. Transferred to reaction chemistry, the 'information' relates to the individual structures being changed. Calculating the difference ΔH of incoming and outgoing paths (ΔH) is feasible, because the network is directed, i.e. links are not bi-directional.
Positive values of Δ describe more 'chemical information' (structures forming the species node) is being transferred to the species than carried away (other species formed via the species node) and vice versa. Figure S5: Visualization of the difference in harmonic centrality of in-and outgoing links.
Unsurprisingly, ehis the species with the lowest value of ΔH, meaning it is involved in many more reactions that consume than yield solvated electrons. O2 and Clare ranked highest, i.e.
fast reactions appear to dominate its formation in comparison to reactions depleting the species' concentration. Interestingly, H2 and H2O2 show a negative ΔH, although high concentrations are predicted by the steady states of the kinetic simulations. The latter are assumed to be attributed to the primary generation rate during radiolysis. This can be vastly underrepresented by the introduced pseudo-reactions, depending on the dose rate.
Moreover, it must be borne in mind that the graph analysis shown here does not take differences in concentration into account but only analyzes the network structure itself. Therefore, species can show node sizes in Figure S4 or ranks in Figure S5 that do not agree with the steady state concentrations simulated. Thus, graph analysis cannot substitute proper kinetic simulation to predict actual chemical conditions during LP-TEM. Instead, it provides an understanding of the underlying connections independent of dose rate and initial concentrations.
A tabular representation of the kinetic model can be found below ( Table S2). The model ensures mass conservation during the temporal evolution in the simulation ( Figure S6).    Figure S8(a) displays a gold nanostructure growth at "sub-threshold" conditions. The process can be partitioned into four phases: First, a rapid agglomeration of (mainly) rod-shaped gold nanoparticles is observed within the first few frames of the experiment. The particles are already present prior to irradiation, which is discussed below. Second, rapid anisotropic growth in the direction of the long axis of the primary particles occurs. The growth rate in terms of projected area per time t in this phase amounts to (100 ± 3) e -(nm²s) -1 (see orange line in Figure S8(b)). This is comparable to the findings obtained by Park et al. [37] at 4000 e -(nm²s) -1 . There, however, a strong shapedependence of the growth rate is reported, which could explain the rapid growth observed here. Third, a substantially slower gold growth is evident (13.1±0.3) e -(nm²s) -1 , green line in Figure S8(b), which is accompanied by a change in the growth direction of about 120°, which fits well to the sixfold fcc symmetry of gold nanoparticles. Growth and reshaping are illustrated in Figure S8(c). Last, after about 52 s, the projected area remains approximately constant, while the crystal structure undergoes slight morphology changes only. The decrease in growth rate suggests that the local concentration gradients are seized so that the amount of solid gold remains in a dynamic equilibrium with the surrounding solution.
Low electron flux density-mediated nucleation has been observed repetitively, as shown below:  S2.2. Additional data on X-ray induced gold precipitation Note that the anisotropic particle distribution is due to the X-ray beam was not irradiating the whole sample area. This is caused by an anisotropic beam-shaping aperture that is installed prior to the liquid cell. Additionally, diffusion could contribute to the asymmetric shape of the area where precipitates are visible.

S4 Origin of the primary gold particles
During our experiments performed in GSMLCs, gold particles have been present already at first sight, even when broadcasting through the wells. This is remarkable, because no particles adhered to the graphene membranes without irradiation, as investigated by control experiments.
The reduction of gold was verified via HRTEM measurement of the lattice plane distances using diffractograms ( Figure S13).
Moreover, we did not observe a comparable experiment in liquid cells comprising two silicon nitride membranes, such as static liquid cells and flow cells. Here, only pattern formation and particle growth in direct vicinity to the irradiated area was observed within minutes, as reported earlier. [92]  Thus, we hypothesize that the presence of the gold nanoparticles is related to the graphene membrane. Negative charging of thin carbon-based membranes during TEM has been reported, in particular of amorphous carbon (aC) [89,90] which can transfer onto GSMLCs during sealing. [69] In addition with the enhanced electrical conductivity of graphene compared to aC, [91] this could lead to a charge dissipation into regions not yet investigated during LP-TEM, spreading possible reduction sites even across different microwells. Moreover, this theory offers a possible explanation towards the mechanisms behind apparent mitigation of radiolysis and beam damage in graphene-based liquid cells. [26,27] To test our theory, we mimicked possible electron beam-induced charging by conducting galvanostatic cathodic polarizations of different currents and durations applied in 20 mм HAuCl4 solution. TEM grids covered with graphene were harnessed as the working electrode.
The grids were prepared identical to those used for GSMLC-sealing. A standard three-electrode system was used with a leakless Ag/AgCl electrode (ET072-1, eDAQ Pty Ltd.) as reference electrode and a Pt sheet as a counter electrode. A Zahner Zennium Electrochemical Workstation was employed for electrochemical tests.
At a low cathodic current of 200 nA, nanoparticles have been forming on the sample already after 5 min (Figure S14(a)). HRTEM investigations (Figure S14(b)) revealed i. a. multiple twinned structures matching well to primary particles etched in situ during LP-TEM ( Figure S13(a) By applying the same cathodic currents to pristine GSMLC carrier frames, reference experiments were performed on silicon nitride windows. Afterwards, no particles were observed via TEM on the free-standing membranes. This is, however, not surprising due to the insulating nature of silicon nitride.   Assuming a spherical particle geometry to estimate the surface of the particle A is only a coarse approximation (c.f. the particle in Figure 6(a)) of the true surface area A'. If shape evolutions are occurring isotropically, the related misfit of the model is correctable with a constant shapetransformation factor s. Hence, by regarding only changes of isotropic nature relative to an initial surface A0 this factor is cancelled out and the spherical approximation remains applicable.

S7 Concentration-dependent nanoparticle dissolution
The reduction of mass m of a nanoparticle over time t can be described via the Noyes-Whitney equation, with respect to the diffusion coefficient of the solvated species D, the surface area of the particle A, the liquid thickness of the boundary layer surrounding the particle L, the solubility of species cs, and the equilibrium concentration of the dissolved species ∞ : [95] = ( − ∞ ) (S5) Assuming a spherical particle with density ρ surrounded by a highly super-saturated solution allows reshaping of (S5) by expressing the mass in terms of density and radius and assuming ∞ being small compared to cs: [31,95] = Equation (S9) can be fitted using a least-squares optimization of the parameters used for solving the ODE. In particular, this denotes the empirical pre-factor B, the boundary condition of r at the beginning of observation (rt=0) and the desired quantity rcrit. This is exemplarily shown in Figure S19 for the primary particle in Figure S12.
Figure S19: Exemplary fit of the critical radius on the particle shown in Figure S12 corresponding to the highest R 2 adj. in Figure S20.
Consequently, the numerical ODE solver is connected in series with the least square optimization algorithm. As both tools converge against (possible) local minima, the goodness of the outcome (quantified by the adjusted coefficient of determination R 2 adj.) is highly sensitive to the starting parameters. Therefore, those were varied systematically to detect the global minimum within this parameter space ( Figure S20). Figure S20: Exemplary parameter space in critical radius fitting. The colorbar illustrates the adjusted coefficient of determination R 2 adj..

S8 Chemical potential calculation
In the following, the transition from a purely solid gold phase into partially dissolved gold in steady state condition with the residual solid phase after infinite time is discussed: Au solid → Au solid + Au solvated (S10) Mass conservation demands that the total amount of Au remains constant: Au solid ( = 0) = Au solid ( = ∞) + Au solvated ( = ∞) (S11) As discussed above, stray irradiation, irradiation during alignment, and graphene-modulated charge-dissipation reduced a major amount of the gold initially stored in HAuCl4 to Au 0 . As the volume remains unchanged, this can be directly expressed by concentrations c, where c0 is the initially solvated amount.
If the transition does reach a steady state situation, the net concentration in the phases does not change anymore. Thus, the chemical potentials of these phases must be equalized. Therefore, describing any of the final phases is sufficient for calculating the chemical potential difference . At a given thermal energy B , of the transition shown in (S10) is characterized by the fraction of solid and solvated amount of gold prior and after the transition, 0 and ∞ . [93,94] = ( 0 ∞ ) (S14) Accounting for the description of the two situations above yields the following expressions for Inserting (S15) and (S16) in (S14) yields the desired expression shown in the main manuscript: ( 1 steady state (Au 0 ) 0 (HAuCl 4 ) ) = B ln ( 0 (HAuCl 4 ) steady state (Au 0 ) ) (S17)

S9 Supporting video captions
Video S1: In situ LP-TEM video showing gold nanoparticle dynamics in 20 mм HAuCl4 under repetitive dose-rate change. Speed is 5-fold accelerated. A scale bar can be found in Figure 6 in the main manuscript.